

rm(list = ls())

###############################################################################
########################### REPLICATION MASTER SCRIPT #########################
###############################################################################

# Article: "Initiate and Elevate! How Political Parties Can Set an Agenda"
# Author: Daniel Sandvej Eriksen
# Replication Master Script

###############################################################################
############################## Set Working Directory ##########################
###############################################################################



# Detect path of this script (works in RStudio, Rscript console, Dataverse)
this_file <- function() {
  cmdArgs <- commandArgs(trailingOnly = FALSE)
  fileFlag <- "--file="
  match <- grep(fileFlag, cmdArgs)
  
  # Running via Rscript from terminal or Dataverse
  if (length(match) > 0) {
    return(normalizePath(sub(fileFlag, "", cmdArgs[match])))
  }
  
  # Running inside RStudio
  if (!is.null(sys.calls())) {
    for (i in sys.calls()) {
      if (all(grepl("source", i, fixed = TRUE)))  {
        return(normalizePath(as.character(attr(i, "srcfile"))))
      }
    }
  }
  
  # Fallback: try RStudio API (if available)
  if ("rstudioapi" %in% installed.packages()[,1]) {
    if (rstudioapi::isAvailable()) {
      p <- rstudioapi::getActiveDocumentContext()$path
      if (!identical(p, "")) return(normalizePath(p))
    }
  }
  
  # If everything fails: stop with a clear error
  stop("Cannot detect script location automatically. Please run with Rscript.")
}

# Set working directory to the folder containing this script
script_path <- this_file()
setwd(dirname(script_path))

cat("Working directory set to:\n", getwd(), "\n")


###############################################################################
############################### Save Session Info #############################
###############################################################################

# Saves full R session information to a file for reproducibility
writeLines(
  capture.output(sessionInfo()),
  file.path(getwd(), "sessionInfo.txt")
)


###############################################################################
############################## Install and load Required Libraries ############
###############################################################################

required_packages <- c(
  "readxl", "openxlsx", "writexl", "dplyr", "ggplot2", "scales", "grid",
  "gridExtra", "patchwork", "tidyr", "pscl", "MASS", "margins", "sjPlot"
)

# Install missing packages
missing_packages <- required_packages[!(required_packages %in% installed.packages()[,"Package"])]

if (length(missing_packages) > 0) {
  cat("Installing missing packages:\n")
  print(missing_packages)
  install.packages(missing_packages, dependencies = TRUE)
}

# Load all packages
for (pkg in required_packages) {
  library(pkg, character.only = TRUE)
}


###############################################################################
############################### Define Paths ##################################
###############################################################################

# Defines the paths to the input data files and output directory.
data_path    <- "../data/analysis"
results_dir  <- "../results"

if (!dir.exists(results_dir)) {
  dir.create(results_dir, recursive = TRUE)
}

###############################################################################
############################### Load Data #####################################
###############################################################################

# Table 2
load(file.path(data_path, "data_table_2.Rdata"))

# Table 3
load(file.path(data_path, "data_table_3.Rdata"))

# Figure 2
load(file.path(data_path, "data_figure_2.Rdata"))

# Figure 3
load(file.path(data_path, "data_figure_3.Rdata"))

# Figure 4
load(file.path(data_path, "data_figure_4.Rdata"))

# Figure 5
load(file.path(data_path, "data_figure_5.Rdata"))

# Figure 6
load(file.path(data_path, "data_figure_6.Rdata"))

# Figure 7
load(file.path(data_path, "data_figure_7.Rdata"))

# Figure 8
load(file.path(data_path, "data_figure_8.Rdata"))

# Figure 9
load(file.path(data_path, "data_figure_9.Rdata"))

# Figure 10
load(file.path(data_path, "data_figure_10.Rdata"))

# Figures 11–12
load(file.path(data_path, "data_figure_11_12.Rdata"))

# Figure 13
load(file.path(data_path, "data_figure_13.Rdata"))

# Figure A1
load(file.path(data_path, "data_figure_a1.Rdata"))

# Figure A2
load(file.path(data_path, "data_figure_a2.Rdata"))

# Figure A3
load(file.path(data_path, "data_figure_a3.Rdata"))

# Figure A4
load(file.path(data_path, "data_figure_a4.Rdata"))

# Figure A5
load(file.path(data_path, "data_figure_a5.Rdata"))

# Figure A6
load(file.path(data_path, "data_figure_a6.Rdata"))

# Figure A7
load(file.path(data_path, "data_figure_a7.Rdata"))

# Figure A8
load(file.path(data_path, "data_figure_a8.Rdata"))

# Figure A9
load(file.path(data_path, "data_figure_a9.Rdata"))

# Figure A10
load(file.path(data_path, "data_figure_a10.Rdata"))

# Figure A11
load(file.path(data_path, "data_figure_a11.Rdata"))

# Figure A12
load(file.path(data_path, "data_figure_a12.Rdata"))

# Figure A13
load(file.path(data_path, "data_figure_a13.Rdata"))

# Figure A14
load(file.path(data_path, "data_figure_a14.Rdata"))

# Figure A15
load(file.path(data_path, "data_figure_a15.Rdata"))

# Figure A16
load(file.path(data_path, "data_figure_a16.Rdata"))

# Figure A17
load(file.path(data_path, "data_figure_a17.Rdata"))

# Figure A18
load(file.path(data_path, "data_figure_a18.Rdata"))

# Figure A19
load(file.path(data_path, "data_figure_a19.Rdata"))

# Figure A20
load(file.path(data_path, "data_figure_a20.Rdata"))

# Figure A21
load(file.path(data_path, "data_figure_a21.Rdata"))

# Figure A22
load(file.path(data_path, "data_figure_a22.Rdata"))

# Figure A23
load(file.path(data_path, "data_figure_a23.Rdata"))

# Figure A24
load(file.path(data_path, "data_figure_a24.Rdata"))

# Figure A25
load(file.path(data_path, "data_figure_a25.Rdata"))

# Figure A26
load(file.path(data_path, "data_figure_a26.Rdata"))

# Models A1–A3
load(file.path(data_path, "data_models_a1_a2_a3.Rdata"))

# Models A4–A7
load(file.path(data_path, "data_models_a4_a5_a6_a7.Rdata"))

# Model A8
load(file.path(data_path, "data_model_a8.Rdata"))

# Models A9–A17
load(file.path(data_path, "data_models_a9_a10_a11_a12_a13_a14_a15_a16_a17.Rdata"))

# Model A18
load(file.path(data_path, "data_model_a18.Rdata"))

# Model A19
load(file.path(data_path, "data_model_a19.Rdata"))

# Models A20–A21
load(file.path(data_path, "data_models_a20_a21.Rdata"))

# Models A22–A23
load(file.path(data_path, "data_model_a22_a23.Rdata"))

# Models A24–A25
load(file.path(data_path, "data_models_a24_a25.Rdata"))

# Models A26-A27
load(file.path(data_path, "data_models_a26_a27.Rdata"))

# Models A28–A30
load(file.path(data_path, "data_models_a28_a29_a30.Rdata"))

# Model A31
load(file.path(data_path, "data_model_a31.Rdata"))

# Models A32–A35
load(file.path(data_path, "data_models_a32_a33_a34_a35.Rdata"))

# Models A36–A37
load(file.path(data_path, "data_models_a36_a37.Rdata"))

# Table A2
load(file.path(data_path, "data_table_a2.Rdata"))

# Share of initiations about same issue on the same day by the same party (p. 15 in the manuscript)
load(file.path(data_path, "data_initiation_same_issue_same_party_same_day.Rdata"))

# Average days between initiations about same the issue by the same party 
load(file.path(data_path, "data_average_time_between_initiations.Rdata"))



######################################################## Manuscript ####################################################################
########################################################################################################################################


###############################################################################
############################### Table 2 #######################################
###############################################################################

#### Output
write.csv(
  data_table_2,
  file = file.path(results_dir, "Table_2.csv"),
  row.names = FALSE
)


###############################################################################
############################### Table 3 #######################################
###############################################################################


#### Calculating
dk_mean <- round(mean(data_table_3$engagement_share[data_table_3$country == "DK"], na.rm = TRUE), 1)
uk_mean <- round(mean(data_table_3$engagement_share[data_table_3$country == "UK"], na.rm = TRUE), 1)
across_countries <- round(mean(data_table_3$engagement_share, na.rm = TRUE), 1)

table_3 <- data.frame(
  Country = c("The UK", "Denmark", "Across countries"),
  `Average share (%)` = c(uk_mean, dk_mean, across_countries),
  check.names = FALSE
)


# Save as Excel
wb <- createWorkbook()
addWorksheet(wb, "Table 3")
writeData(wb, sheet = "Table 3", x = table_3, startCol = 1, startRow = 1)
saveWorkbook(wb,
             file = file.path(results_dir, "Table_3.xlsx"),
             overwrite = TRUE)



###############################################################################
############################### Figure 2 ######################################
###############################################################################

#### UK subset
figure_2_uk_temp <- ggplot(
  subset(data_figure_2, country == "UK"),
  aes(x = date, y = initiation_count, fill = party)
) +
  geom_bar(stat = "identity") +
  scale_fill_manual(values = c("Labour" = "red", "Conservatives" = "blue"),
                    na.value = "gray") +
  theme_minimal() +
  theme(axis.text.x = element_text(size = 20, angle = 90),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 22),
        axis.title.y = element_text(size = 22),
        strip.text = element_text(size = 22),
        legend.position = "none",
        plot.title = element_text(size = 24, hjust = 0.5)) +
  labs(x = "", y = "Frequency\n", title = "UK") +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 8),
                     limits = c(0, 14))

#### DK subset
figure_2_dk_temp <- ggplot(
  subset(data_figure_2, country == "DK"),
  aes(x = date, y = initiation_count, fill = party)
) +
  geom_bar(stat = "identity") +
  scale_fill_manual(values = c("Socialdemokratiet" = "red", "Venstre" = "blue"),
                    na.value = "gray") +
  theme_minimal() +
  theme(axis.text.x = element_text(size = 20, angle = 90),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 22),
        axis.title.y = element_text(size = 22),
        strip.text = element_text(size = 22),
        legend.position = "none",
        plot.title = element_text(size = 24, hjust = 0.5)) +
  labs(x = "\nDate", y = "Frequency\n", title = "DK") +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 8),
                     limits = c(0, 14))


#### Combine vertically
figure_2 <- grid.arrange(
  figure_2_uk_temp,
  figure_2_dk_temp,
  ncol = 1
)

#### Output
ggsave(
  filename = "figure_2.pdf",
  plot = figure_2,
  path = results_dir,
  width = 10,
  height = 12
)

###############################################################################
############################### Figure 3 ######################################
###############################################################################

# UK
figure_3_uk_temp <- ggplot(
  subset(data_figure_3, country == "UK"),
  aes(x = days_before,
      y = average_tweet_count_all_other_parties,
      group = group,
      linetype = group)
) +
  geom_line(size = 1.5, color = "black") +
  geom_point(size = 3, color = "black") +
  labs(
    title = "UK",
    x = "\nDays before the tweet",
    y = "Average number of tweets\n"
  ) +
  theme_minimal() +
  ylim(0, 150) +
  scale_linetype_manual(values = c("Initiations" = "solid",
                                   "Non-Initiations" = "dashed")) +
  theme(
    plot.title = element_text(size = 18, hjust = 0.5),
    text = element_text(size = 18),
    legend.position = "none",
    axis.text.x = element_text(size = 18, margin = margin(t = 20), angle = 25),
    axis.text.y = element_text(size = 18, margin = margin(r = 20)),
    axis.title  = element_text(size = 20, margin = margin(b = 20))
  )


# DK
figure_3_dk_temp <- ggplot(
  subset(data_figure_3, country == "DK"),
  aes(x = days_before,
      y = average_tweet_count_all_other_parties,
      group = group,
      linetype = group)
) +
  geom_line(size = 1.5, color = "black") +
  geom_point(size = 3, color = "black") +
  labs(
    title = "DK",
    x = "\nDays before the tweet",
    y = "Average number of tweets\n"
  ) +
  theme_minimal() +
  ylim(0, 20) +
  scale_linetype_manual(values = c("Initiations" = "solid",
                                   "Non-Initiations" = "dashed")) +
  theme(
    plot.title = element_text(size = 18, hjust = 0.5),
    text = element_text(size = 18),
    legend.position = "none",
    axis.text.x = element_text(size = 18, margin = margin(t = 20), angle = 25),
    axis.text.y = element_text(size = 18, margin = margin(r = 20)),
    axis.title  = element_text(size = 20, margin = margin(b = 20))
  )


# Combine 
figure_3 <- wrap_plots(
  figure_3_uk_temp, 
  figure_3_dk_temp,
  ncol = 2,
  widths = c(1, 1)
)

ggsave(
  filename = "figure_3.pdf",
  plot = figure_3,
  path = results_dir,
  width = 15,
  height = 6
)

###############################################################################
############################### Figure 4 ######################################
###############################################################################


# UK
figure_4_uk_temp <- ggplot(
  subset(data_figure_4, country == "UK"),
  aes(x = days_before,
      y = mean_value,
      group = group,
      linetype = group)
) +
  geom_line(size = 1.5, color = "black") +
  geom_point(size = 3, color = "black") +
  labs(
    title = "UK",
    x = "\nDays before the tweet",
    y = "Average number of news articles\n"
  ) +
  theme_minimal() +
  ylim(0, 30) +
  scale_linetype_manual(values = c("Initiations" = "solid",
                                   "Non-Initiations" = "dashed")) +
  theme(
    plot.title = element_text(size = 18, hjust = 0.5),
    text = element_text(size = 18),
    legend.position = "none",
    axis.text.x = element_text(size = 18, margin = margin(t = 20), angle = 25),
    axis.text.y = element_text(size = 18, margin = margin(r = 20)),
    axis.title  = element_text(size = 20, margin = margin(b = 20))
  )


# DK
figure_4_dk_temp <- ggplot(
  subset(data_figure_4, country == "DK"),
  aes(x = days_before,
      y = mean_value,
      group = group,
      linetype = group)
) +
  geom_line(size = 1.5, color = "black") +
  geom_point(size = 3, color = "black") +
  labs(
    title = "DK",
    x = "\nDays before the tweet",
    y = "Average number of news articles\n"
  ) +
  theme_minimal() +
  ylim(0, 10) +
  scale_linetype_manual(values = c("Initiations" = "solid",
                                   "Non-Initiations" = "dashed")) +
  theme(
    plot.title = element_text(size = 18, hjust = 0.5),
    text = element_text(size = 18),
    legend.position = "none",
    axis.text.x = element_text(size = 18, margin = margin(t = 20), angle = 25),
    axis.text.y = element_text(size = 18, margin = margin(r = 20)),
    axis.title  = element_text(size = 20, margin = margin(b = 20))
  )


# Combine
figure_4 <- wrap_plots(
  figure_4_uk_temp, 
  figure_4_dk_temp,
  ncol = 2,
  widths = c(1, 1)
)


ggsave(
  filename = "figure_4.pdf",
  plot = figure_4,
  path = results_dir,
  width = 15,
  height = 8
)

###############################################################################
############################### Figure 5 ######################################
###############################################################################


# UK
figure_5_uk_temp <- ggplot(
  subset(data_figure_5, country == "UK"),
  aes(x = day,
      y = average_elevation_tweets_count,
      group = group,
      linetype = group)
) +
  geom_line(size = 1.0, color = "black") +
  geom_point(size = 3, color = "black") +
  labs(
    title = "UK",
    x = "\nDay",
    y = "Average number of tweets about the issue of the post\n"
  ) +
  theme_minimal() +
  scale_linetype_manual(values = c("Non-initiation" = "dashed",
                                   "Initiation"     = "solid")) +
  scale_x_discrete(
    limits = rev(levels(subset(data_figure_5, country == "UK")$day))
  ) +
  theme(
    text = element_text(size = 18),
    legend.position = "none",
    axis.text.x = element_text(size = 18, margin = margin(t = 22)),
    axis.text.y = element_text(size = 18, margin = margin(r = 22)),
    axis.title  = element_text(size = 20, margin = margin(b = 22))
  ) +
  geom_text(
    aes(label = round(average_elevation_tweets_count, 1)),
    vjust = -0.5,
    size = 6
  )


# DK
figure_5_dk_temp <- ggplot(
  subset(data_figure_5, country == "DK"),
  aes(x = day,
      y = average_elevation_tweets_count,
      group = group,
      linetype = group)
) +
  geom_line(size = 1.0, color = "black") +
  geom_point(size = 3, color = "black") +
  labs(
    title = "DK",
    x = "\nDay",
    y = "Average number of tweets about the issue of the post\n"
  ) +
  theme_minimal() +
  scale_linetype_manual(values = c("Non-initiation" = "dashed",
                                   "Initiation"     = "solid")) +
  scale_x_discrete(
    limits = rev(levels(subset(data_figure_5, country == "DK")$day))
  ) +
  theme(
    text = element_text(size = 18),
    legend.position = "none",
    axis.text.x = element_text(size = 18, margin = margin(t = 22)),
    axis.text.y = element_text(size = 18, margin = margin(r = 22)),
    axis.title  = element_text(size = 20, margin = margin(b = 22))
  ) +
  geom_text(
    aes(label = round(average_elevation_tweets_count, 1)),
    vjust = -0.5,
    size = 6
  )


# Combine 
figure_5 <- wrap_plots(
  figure_5_uk_temp, 
  figure_5_dk_temp,
  ncol = 2,
  widths = c(1, 1)
)


ggsave(
  filename = "figure_5.pdf",
  plot = figure_5,
  path = results_dir,
  width = 12,
  height = 8
)



###############################################################################
############################### Figure 6 ######################################
###############################################################################

#### Plotting
figure_6 <- ggplot(data_figure_6,
                   aes(x = elevation_count_on_initiation_day, y = n)) +
  geom_bar(stat = "identity") +
  geom_vline(xintercept = data_figure_6$mean_elevation,
             linetype = "dotted", color = "green", size = 2) +
  theme_minimal() +
  theme(axis.text.x = element_text(size = 20, angle = 90),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 22),
        axis.title.y = element_text(size = 22),
        strip.text = element_text(size = 22),
        legend.position = "none") +
  labs(x = "\nElevation ", y = "Frequency\n")

#### Output figure
ggsave(
  filename = "figure_6.pdf",
  plot = figure_6,
  path = results_dir,
  width = 10,
  height = 8
)


#### Calculate and output mean and standard deviation
descriptives_figure_6 <- data.frame(
  mean_elevation = round(unique(data_figure_6$mean_elevation), 1),
  sd_elevation   = round(unique(data_figure_6$sd_elevation), 1)
)
write.csv(
  descriptives_figure_6,
  file = file.path(results_dir, "Figure_6_descriptives.csv"),
  row.names = FALSE
)


###############################################################################
############################### Figure 7 ######################################
###############################################################################

#### Plot
figure_7 <- data_figure_7 %>%
  ggplot(aes(x = period, y = irf_elevation_reaction,
             ymin = lower_elevation_reaction, ymax = upper_elevation_reaction)) +
  geom_hline(yintercept = 0, color = "red") +
  geom_ribbon(fill = "grey", alpha = 0.2) +
  geom_line() +
  theme_light() +
  ylab("Number of tweets by competing party MPs\n") +
  xlab("\nNumber of minutes since the one-tweet increase in elevation") +
  theme_minimal() +
  theme(axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 22),
        axis.title.y = element_text(size = 22),
        strip.text = element_text(size = 22)) +
  scale_x_continuous(breaks = data_figure_7$period)

#### Output
ggsave(
  filename = "figure_7.pdf",
  plot = figure_7,
  path = results_dir,
  width = 12,
  height = 8
)


###############################################################################
############################### Figure 8 ######################################
###############################################################################

#### Plot
figure_8 <- data_figure_8 %>%
  ggplot(aes(x = period, y = irf_elevation_reaction,
             ymin = lower_elevation_reaction, ymax = upper_elevation_reaction)) +
  geom_hline(yintercept = 0, color = "red") +
  geom_ribbon(fill = "grey", alpha = 0.2) +
  geom_line() +
  theme_light() +
  ylab("Cumulated number of tweets by competing party MPs\n") +
  xlab("\nNumber of minutes since the one-tweet increase in elevation") +
  theme_minimal() +
  theme(axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 18.2),
        axis.title.y = element_text(size = 18.2),
        strip.text = element_text(size = 22)) +
  scale_x_continuous(breaks = data_figure_8$period)

#### Output
ggsave(
  filename = "figure_8.pdf",
  plot = figure_8,
  path = results_dir,
  width = 12,
  height = 9
)


###############################################################################
############################### Figure 9 ######################################
###############################################################################

#### Run model
model_for_figure_9_temp <- glm.nb(
  formula = engagement_tweets_following_week ~ elevation_tweets +
    issue_salience_last_week_average +
    factor(year) +
    factor(country) +
    propose_solution,
  data = data_figure_9
)

#### Output
output_path <- file.path(results_dir, "figure_9.pdf")

pdf(file = output_path, width = 10, height = 8)

par(mar = c(4, 6.5, 2, 2))

figure_9 <- cplot(
  model_for_figure_9_temp,
  "elevation_tweets",
  draw  = TRUE,
  main  = "",
  xlab  = "Elevation",
  ylab  = "",
  cex.axis = 1.4,
  cex.lab  = 1.6,
  cex.main = 1.6,
  cex.sub  = 1.3
)

mtext("Predicted number of engagement tweets in the following week",
      side = 2, line = 5, cex = 1.6)

dev.off()

###############################################################################
############################### Figure 10 #####################################
###############################################################################

#### Run model
model_for_figure_10_temp <- glm.nb(
  formula = news_articles_day_after ~ elevation_tweets +
    issue_salience_last_week_average +
    factor(year) +
    factor(country) +
    propose_solution,
  data = data_figure_10
)

#### Output
output_path <- file.path(results_dir, "figure_10.pdf")

pdf(file = output_path, width = 10, height = 8)

figure_10 <- cplot(
  model_for_figure_10_temp, "elevation_tweets", draw = TRUE,
  main = "",
  xlab = "Elevation",
  ylab = "Predicted number of news articles on the following day",
  cex.axis = 1.4,
  cex.lab  = 1.5,
  cex.main = 1.7,
  cex.sub  = 1.2
)

dev.off()


###############################################################################
############################### Figure 11 #####################################
###############################################################################

#### Run model
model_for_figure_11_temp <- glm.nb(
  formula = news_articles_week_after ~ elevation_tweets +
    issue_salience_last_week_average +
    factor(year) +
    factor(country) +
    propose_solution,
  data = data_figures_11_12
)

#### Output
output_path <- file.path(results_dir, "figure_11.pdf")

pdf(file = output_path, width = 10, height = 8)

figure_11 <- cplot(
  model_for_figure_11_temp, "elevation_tweets", draw = TRUE,
  main = "",
  xlab = "Elevation",
  ylab = "Predicted number of news articles in the following week",
  cex.axis = 1.4,
  cex.lab  = 1.5,
  cex.main = 1.7,
  cex.sub  = 1.2
)

dev.off()


###############################################################################
############################### Figure 12 #####################################
###############################################################################

#### Run model
model_for_figure_12_temp <- glm.nb(
  formula = news_articles_week_after ~ engagement_tweets +
    issue_salience_last_week_average +
    factor(year) +
    factor(country) +
    propose_solution,
  data = data_figures_11_12
)

#### Output
output_path <- file.path(results_dir, "figure_12.pdf")

pdf(file = output_path, width = 10, height = 8)

figure_12 <- cplot(
  model_for_figure_12_temp, "engagement_tweets", draw = TRUE,
  main = "",
  xlab = "Engagement",
  ylab = "Predicted number of news articles in the following week",
  cex.axis = 1.4,
  cex.lab  = 1.5,
  cex.main = 1.7,
  cex.sub  = 1.2
)

dev.off()


###############################################################################
############################### Figure 13 #####################################
###############################################################################

#### Run model
model_for_figure_13_temp <- zeroinfl(
  total_followup_questions ~ elevation_tweets +
    factor(country) +
    factor(year),
  data = data_figure_13,
  dist = "negbin"
)

#### Output
count_coef_figure_13 <- summary(model_for_figure_13_temp)$coefficients$count["elevation_tweets", ]

irr_data_figure_13 <- data.frame(
  Component  = "",
  Estimate   = exp(count_coef_figure_13["Estimate"]),
  SE         = count_coef_figure_13["Std. Error"],
  Lower      = exp(count_coef_figure_13["Estimate"] - 1.96 * count_coef_figure_13["Std. Error"]),
  Upper      = exp(count_coef_figure_13["Estimate"] + 1.96 * count_coef_figure_13["Std. Error"]),
  Significant = count_coef_figure_13["Pr(>|z|)"] < 0.05
)

figure_13 <- ggplot(
  irr_data_figure_13,
  aes(x = Component, y = Estimate, color = Significant, fill = Significant)
) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "gray50") +
  geom_errorbar(aes(ymin = Lower, ymax = Upper), width = 0.2, size = 1) +
  geom_point(size = 4, shape = 21, stroke = 1.5) +
  scale_color_manual(values = c("FALSE" = "gray60", "TRUE" = "black")) +
  scale_fill_manual(values  = c("FALSE" = "white",  "TRUE" = "black")) +
  theme_minimal() +
  theme(
    legend.position = "none",
    axis.text.x = element_text(size = 14),
    axis.text.y = element_text(size = 14),
    axis.title.x = element_text(size = 18),
    axis.title.y = element_text(size = 18)
  ) +
  labs(
    title    = "",
    subtitle = "",
    x        = "Elevation",
    y        = "Incidence Rate Ratio"
  )

ggsave(
  filename = "figure_13.pdf",
  plot = figure_13,
  path = results_dir,
  width = 12,
  height = 8
)



###############################################################################
############################### Page 15 (manuscript), 6.2% ####################
###############################################################################

#### Calculate
duplicate_combinations <- data_initiation_same_issue_same_party_same_day %>%
  group_by(party, CAP_issue, day) %>%
  summarise(count = n(), .groups = 'drop') %>%
  filter(count > 1) %>%
  arrange(desc(count))

value_6.2_page_15 <- sum(duplicate_combinations$count) / nrow(data_initiation_same_issue_same_party_same_day) * 100
value_6.2_page_15 <- round(value_6.2_page_15, 1)

#### Output
initiation_same_issue_same_party_same_day_value <- data.frame(
  `Share (%)` = value_6.2_page_15,
  check.names = FALSE
)

write.csv(
  initiation_same_issue_same_party_same_day_value,
  file = file.path(results_dir, "Initiation_same_issue_same_party_same_day.csv"),
  row.names = FALSE
)


###############################################################################
############################### Page 15 (manuscript), 100 days ################
###############################################################################

#### Calculate
average_time_gaps <- data_average_time_between_initiations %>%
  arrange(party, CAP_issue, day) %>%
  group_by(party, CAP_issue) %>%
  mutate(
    previous_day = lag(day),
    time_gap_days = as.numeric(day - previous_day)
  ) %>%
  filter(!is.na(time_gap_days))

average_time_gaps <- average_time_gaps[average_time_gaps$time_gap_days > 0,]
average_time_gaps_value <- mean(average_time_gaps$time_gap_days)
average_time_gaps_value <- round(average_time_gaps_value, 1)



#### Output
average_time_gaps_between_initiations <- data.frame(
  `Days` = average_time_gaps_value,
  check.names = FALSE
)

write.csv(
  average_time_gaps_between_initiations,
  file = file.path(results_dir, "Average_days_between_initiations_about_same_issue_by_same_party.csv"),
  row.names = FALSE
)






#################################################### Appendices ##################################################################
##################################################################################################################################


###############################################################################
############################### Model A1 ######################################
###############################################################################

#### Producing
model_a1 <- glm.nb(
  formula = engagement_tweets_following_week ~ elevation_tweets +
    issue_salience_last_week_average +
    factor(year) +
    factor(country) +
    propose_solution,
  data = data_models_a1_a2_a3
)

#### Outputting
tab_model(
  model_a1,
  vcov.fun = "HC3" ,
  digits = 4,
  file   = file.path(results_dir, "model_A1.doc")
)


###############################################################################
############################### Model A2 ######################################
###############################################################################

model_a2 <- glm.nb(
  formula = engagement_tweets_following_week ~ elevation_tweets +
    issue_salience_last_week_average +
    factor(year) +
    factor(party) +
    propose_solution,
  data = subset(data_models_a1_a2_a3, country == "uk")
)

#### Outputting
tab_model(
  model_a2,
  vcov.fun = "HC3" ,
  digits = 4,
  file   = file.path(results_dir, "model_A2.doc")
)


###############################################################################
############################### Model A3 ######################################
###############################################################################

model_a3 <- glm.nb(
  formula = engagement_tweets_following_week ~ elevation_tweets +
    issue_salience_last_week_average +
    factor(year) +
    factor(party) +
    propose_solution,
  data = subset(data_models_a1_a2_a3, country == "dk")
)

#### Outputting
tab_model(
  model_a3,
  vcov.fun = "HC3" ,
  digits = 4,
  file   = file.path(results_dir, "model_A3.doc")
)


###############################################################################
############################### Model A4 ######################################
###############################################################################

#### Run model
model_a4 <- glm.nb(
  formula = news_articles_day_after ~ elevation_tweets +
    issue_salience_last_week_average +
    factor(year) +
    factor(country) +
    propose_solution,
  data = data_models_a4_a5_a6_a7
)

#### Outputting
tab_model(
  model_a4,
  vcov.fun = "HC3" ,
  digits = 4,
  file   = file.path(results_dir, "model_A4.doc")
)


###############################################################################
############################### Model A5 ######################################
###############################################################################

#### Run model
model_a5 <- glm.nb(
  formula = news_articles_day_after ~ elevation_tweets +
    issue_salience_last_week_average +
    factor(year) +
    factor(party) +
    propose_solution,
  data = subset(data_models_a4_a5_a6_a7, country == "dk")
)

#### Outputting
tab_model(
  model_a5,
  vcov.fun = "HC3" ,
  digits = 4,
  file   = file.path(results_dir, "model_A5.doc")
)


###############################################################################
############################### Model A6 ######################################
###############################################################################

#### Run model
model_a6 <- glm.nb(
  formula = news_articles_day_after ~ elevation_tweets +
    issue_salience_last_week_average +
    factor(year) +
    factor(party) +
    propose_solution,
  data = subset(data_models_a4_a5_a6_a7, country == "uk")
)

#### Outputting
tab_model(
  model_a6,
  vcov.fun = "HC3" ,
  digits = 4,
  file   = file.path(results_dir, "model_A6.doc")
)


###############################################################################
############################### Model A7 ######################################
###############################################################################

#### Run model
model_a7 <- glm.nb(
  formula = news_articles_day_after ~ elevation_tweets * country +
    issue_salience_last_week_average +
    factor(year) +
    propose_solution,
  data = data_models_a4_a5_a6_a7
)

#### Outputting
tab_model(
  model_a7,
  vcov.fun = "HC3" ,
  digits = 4,
  file   = file.path(results_dir, "model_A7.doc")
)


###############################################################################
############################### Model A8 ######################################
###############################################################################

#### Run model
model_a8 <- glm.nb(
  formula = news_articles_day_after ~ elevation_tweets * government_type +
    issue_salience_last_week_average +
    factor(year) +
    propose_solution,
  data = data_model_a8
)

#### Outputting
tab_model(
  model_a8,
  vcov.fun = "HC3" ,
  digits = 4,
  file   = file.path(results_dir, "model_A8.doc")
)


###############################################################################
############################### Model A9 ######################################
###############################################################################

#### Run model
model_a9 <- glm.nb(
  formula = news_articles_week_after ~ elevation_tweets +
    issue_salience_last_week_average +
    factor(year) +
    factor(country) +
    propose_solution,
  data = data_models_a9_a10_a11_a12_a13_a14_a15_a16_a17
)

#### Outputting
tab_model(
  model_a9,
  vcov.fun = "HC3" ,
  digits = 4,
  file   = file.path(results_dir, "model_A9.doc")
)


###############################################################################
############################### Model A10 #####################################
###############################################################################

#### Run model
model_a10 <- glm.nb(
  formula = news_articles_week_after ~ elevation_tweets +
    issue_salience_last_week_average +
    factor(year) +
    factor(party) +
    propose_solution,
  data = subset(data_models_a9_a10_a11_a12_a13_a14_a15_a16_a17, country == "uk")
)

#### Outputting
tab_model(
  model_a10,
  vcov.fun = "HC3" ,
  digits = 4,
  file   = file.path(results_dir, "model_A10.doc")
)


###############################################################################
############################### Model A11 #####################################
###############################################################################

#### Run model
model_a11 <- glm.nb(
  formula = news_articles_week_after ~ elevation_tweets +
    issue_salience_last_week_average +
    factor(year) +
    factor(party) +
    propose_solution,
  data = subset(data_models_a9_a10_a11_a12_a13_a14_a15_a16_a17, country == "dk")
)

#### Outputting
tab_model(
  model_a11,
  vcov.fun = "HC3" ,
  digits = 4,
  file   = file.path(results_dir, "model_A11.doc")
)


###############################################################################
############################### Model A12 #####################################
###############################################################################

#### Run model
model_a12 <- glm.nb(
  formula = news_articles_day_after ~ engagement_tweets +
    issue_salience_last_week_average +
    factor(year) +
    factor(country) +
    propose_solution,
  data = data_models_a9_a10_a11_a12_a13_a14_a15_a16_a17
)

#### Outputting
tab_model(
  model_a12,
  vcov.fun = "HC3" ,
  digits = 4,
  file   = file.path(results_dir, "model_A12.doc")
)


###############################################################################
############################### Model A13 #####################################
###############################################################################

#### Run model
model_a13 <- glm.nb(
  formula = news_articles_day_after ~ engagement_tweets +
    issue_salience_last_week_average +
    factor(year) +
    factor(party) +
    propose_solution,
  data = subset(data_models_a9_a10_a11_a12_a13_a14_a15_a16_a17, country == "uk")
)

#### Outputting
tab_model(
  model_a13,
  vcov.fun = "HC3" ,
  digits = 4,
  file   = file.path(results_dir, "model_A13.doc")
)


###############################################################################
############################### Model A14 #####################################
###############################################################################

#### Run model
model_a14 <- glm.nb(
  formula = news_articles_day_after ~ engagement_tweets +
    issue_salience_last_week_average +
    factor(year) +
    factor(party) +
    propose_solution,
  data = subset(data_models_a9_a10_a11_a12_a13_a14_a15_a16_a17, country == "dk")
)

#### Outputting
tab_model(
  model_a14,
  vcov.fun = "HC3" ,
  digits = 4,
  file   = file.path(results_dir, "model_A14.doc")
)


###############################################################################
############################### Model A15 #####################################
###############################################################################

#### Run model
model_a15 <- glm.nb(
  formula = news_articles_week_after ~ engagement_tweets +
    issue_salience_last_week_average +
    factor(year) +
    factor(country) +
    propose_solution,
  data = data_models_a9_a10_a11_a12_a13_a14_a15_a16_a17
)

#### Outputting
tab_model(
  model_a15,
  vcov.fun = "HC3" ,
  digits = 4,
  file   = file.path(results_dir, "model_A15.doc")
)


###############################################################################
############################### Model A16 #####################################
###############################################################################

#### Run model
model_a16 <- glm.nb(
  formula = news_articles_week_after ~ engagement_tweets +
    issue_salience_last_week_average +
    factor(year) +
    factor(party) +
    propose_solution,
  data = subset(data_models_a9_a10_a11_a12_a13_a14_a15_a16_a17, country == "dk")
)

#### Outputting
tab_model(
  model_a16,
  vcov.fun = "HC3" ,
  digits = 4,
  file   = file.path(results_dir, "model_A16.doc")
)


###############################################################################
############################### Model A17 #####################################
###############################################################################

#### Run model
model_a17 <- glm.nb(
  formula = news_articles_week_after ~ engagement_tweets +
    issue_salience_last_week_average +
    factor(year) +
    factor(party) +
    propose_solution,
  data = subset(data_models_a9_a10_a11_a12_a13_a14_a15_a16_a17, country == "uk")
)

#### Outputting
tab_model(
  model_a17,
  vcov.fun = "HC3" ,
  digits = 4,
  file   = file.path(results_dir, "model_A17.doc")
)


###############################################################################
############################### Model A18 #####################################
###############################################################################

model_a18 <- glm.nb(
  formula = news_articles_day_after ~ elevation_tweets +
    issue_salience_last_week_average +
    factor(year) +
    propose_solution +
    factor(country),
  data = data_model_a18
)

#### Outputting
tab_model(
  model_a18,
  vcov.fun = "HC3" ,
  digits = 4,
  file   = file.path(results_dir, "model_A18.doc")
)


###############################################################################
############################### Model A19 #####################################
###############################################################################

model_a19 <- glm.nb(
  formula = news_articles_week_after ~ elevation_tweets +
    issue_salience_last_week_average +
    factor(year) +
    propose_solution +
    factor(country),
  data = data_model_a19
)

#### Outputting
tab_model(
  model_a19,
  vcov.fun = "HC3" ,
  digits = 4,
  file   = file.path(results_dir, "model_A19.doc")
)


###############################################################################
############################### Model A20 #####################################
###############################################################################

model_a20 <- glm.nb(
  formula = engagement_tweets_day_after ~ elevation_tweets +
    issue_salience_last_week_average +
    factor(year) +
    factor(country) +
    propose_solution,
  data = data_models_a20_a21
)

#### Outputting
tab_model(
  model_a20,
  vcov.fun = "HC3" ,
  digits = 4,
  file   = file.path(results_dir, "model_A20.doc")
)


###############################################################################
############################### Model A21 #####################################
###############################################################################

model_a21 <- glm.nb(
  formula = engagement_tweets_following_week ~ elevation_tweets +
    issue_salience_last_week_average +
    factor(year) +
    factor(country) +
    propose_solution,
  data = data_models_a20_a21
)

#### Outputting
tab_model(
  model_a21,
  vcov.fun = "HC3" ,
  digits = 4,
  file   = file.path(results_dir, "model_A21.doc")
)


###############################################################################
############################### Model A22 #####################################
###############################################################################

model_a22 <- glm.nb(
  formula = news_articles_day_after ~ elevation_tweets +
    issue_salience_last_week_average +
    factor(year) +
    factor(country) +
    propose_solution,
  data = data_model_a22_a23
)

#### Outputting
tab_model(
  model_a22,
  vcov.fun = "HC3" ,
  digits = 4,
  file   = file.path(results_dir, "model_A22.doc")
)


###############################################################################
############################### Model A23 #####################################
###############################################################################

model_a23 <- glm.nb(
  formula = news_articles_week_after ~ elevation_tweets +
    issue_salience_last_week_average +
    factor(year) +
    factor(country) +
    propose_solution,
  data = data_model_a22_a23
)

#### Outputting
tab_model(
  model_a23,
  vcov.fun = "HC3" ,
  digits = 4,
  file   = file.path(results_dir, "model_A23.doc")
)


###############################################################################
############################### Model A24 #####################################
###############################################################################

model_a24 <- glm.nb(
  formula = engagement_tweets_day_after ~ elevation_tweets +
    issue_salience_last_week_average +
    factor(year) +
    factor(country) +
    propose_solution,
  data = data_models_a24_a25
)

#### Outputting
tab_model(
  model_a24,
  vcov.fun = "HC3" ,
  digits = 4,
  file   = file.path(results_dir, "model_A24.doc")
)


###############################################################################
############################### Model A25 #####################################
###############################################################################

model_a25 <- glm.nb(
  formula = engagement_tweets_following_week ~ elevation_tweets +
    issue_salience_last_week_average +
    factor(year) +
    factor(country) +
    propose_solution,
  data = data_models_a24_a25
)

#### Outputting
tab_model(
  model_a25,
  vcov.fun = "HC3" ,
  digits = 4,
  file   = file.path(results_dir, "model_A25.doc")
)


###############################################################################
############################### Model A26 #####################################
###############################################################################

model_a26 <- glm.nb(
  formula = news_articles_day_after ~ elevation_tweets * initiation +
    issue_salience_last_week_average +
    factor(year) +
    propose_solution +
    factor(country),
  data = data_models_a26_a27
)

#### Outputting
tab_model(
  model_a26,
  vcov.fun = "HC3" ,
  digits = 4,
  file   = file.path(results_dir, "model_A26.doc")
)


###############################################################################
############################### Model A27 #####################################
###############################################################################

model_a27 <- glm.nb(
  formula = news_articles_week_after ~ elevation_tweets * initiation +
    issue_salience_last_week_average +
    factor(year) +
    propose_solution +
    factor(country),
  data = data_models_a26_a27
)

#### Outputting
tab_model(
  model_a27,
  vcov.fun = "HC3" ,
  digits = 4,
  file   = file.path(results_dir, "model_A27.doc")
)


###############################################################################
############################### Model A28 #####################################
###############################################################################

#### Run
model_a28 <- zeroinfl(
  total_followup_questions ~ elevation_tweets +
    factor(country) +
    factor(year),
  data = data_models_a28_a29_a30,
  dist = "negbin"
)

#### Output
tab_model(
  model_a28,
  vcov.fun = function(x) sandwich::vcovHC(x, type = "HC3"),
  digits = 4,
  file   = file.path(results_dir, "model_A28.doc")
)



###############################################################################
############################### Model A29 #####################################
###############################################################################

#### Run
model_a29 <- zeroinfl(
  total_followup_questions ~ elevation_tweets +
    factor(year),
  data = subset(data_models_a28_a29_a30, country == "uk"),
  dist = "negbin"
)

#### Output
tab_model(
  model_a29,
  vcov.fun = function(x) sandwich::vcovHC(x, type = "HC3"),
  digits = 4,
  file   = file.path(results_dir, "model_A29.doc")
)



###############################################################################
############################### Model A30 #####################################
###############################################################################

#### Run
model_a30 <- zeroinfl(
  total_followup_questions ~ elevation_tweets +
    factor(year),
  data = subset(data_models_a28_a29_a30, country == "dk"),
  dist = "negbin"
)

#### Output
tab_model(
  model_a30,
  vcov.fun = function(x) sandwich::vcovHC(x, type = "HC3"),
  digits = 4,
  file   = file.path(results_dir, "model_A30.doc")
)



###############################################################################
############################### Model A31 #####################################
###############################################################################

#### Run
model_a31 <- zeroinfl(
  total_later_questions ~ total_early_elevation,
  data = subset(data_model_a31, country == "dk"),
  dist = "negbin"
)

#### Output
tab_model(
  model_a31,
  vcov.fun = function(x) sandwich::vcovHC(x, type = "HC3"),
  digits = 4,
  file   = file.path(results_dir, "model_A31.doc")
)



###############################################################################
############################### Model A32 #####################################
###############################################################################

#### Run
model_a32 <- glm.nb(
  formula = competing_parties_following_day ~ elevation_tweets +
    news_media_salience_last_week_average +
    factor(party) +
    factor(year),
  data = data_models_a32_a33_a34_a35
)

#### Output
tab_model(
  model_a32,
  vcov.fun = "HC3" ,
  digits = 4,
  file   = file.path(results_dir, "model_A32.doc")
)


###############################################################################
############################### Model A33 #####################################
###############################################################################

#### Run
model_a33 <- glm.nb(
  formula = competing_parties_following_week ~ elevation_tweets +
    news_media_salience_last_week_average +
    factor(party) +
    factor(year),
  data = data_models_a32_a33_a34_a35
)

#### Output
tab_model(
  model_a33,
  vcov.fun = "HC3" ,
  digits = 4,
  file   = file.path(results_dir, "model_A33.doc")
)


###############################################################################
############################### Model A34 #####################################
###############################################################################

#### Run
model_a34 <- glm.nb(
  formula = same_bloc_parties_following_day ~ elevation_tweets +
    news_media_salience_last_week_average +
    factor(party) +
    factor(year),
  data = data_models_a32_a33_a34_a35
)

#### Output
tab_model(
  model_a34,
  vcov.fun = "HC3" ,
  digits = 4,
  file   = file.path(results_dir, "model_A34.doc")
)


###############################################################################
############################### Model A35 #####################################
###############################################################################

#### Run
model_a35 <- glm.nb(
  formula = same_bloc_parties_following_week ~ elevation_tweets +
    news_media_salience_last_week_average +
    factor(party) +
    factor(year),
  data = data_models_a32_a33_a34_a35
)

#### Output
tab_model(
  model_a35,
  vcov.fun = "HC3" ,
  digits = 4,
  file   = file.path(results_dir, "model_A35.doc")
)


###############################################################################
############################### Model A36 #####################################
###############################################################################

#### Run
model_a36 <- glm.nb(
  formula = news_articles_day_after ~ elevation_tweets +
    news_articles_last_week_average +
    factor(party) +
    factor(year),
  data = data_models_a36_a37
)

#### Output
tab_model(
  model_a36,
  vcov.fun = "HC3" ,
  digits = 4,
  file   = file.path(results_dir, "model_A36.doc")
)


###############################################################################
############################### Model A37 #####################################
###############################################################################

#### Run
model_a37 <- glm.nb(
  formula = news_articles_week_after ~ elevation_tweets +
    news_articles_last_week_average +
    factor(party) +
    factor(year),
  data = data_models_a36_a37
)

#### Output
tab_model(
  model_a37,
  vcov.fun = "HC3" ,
  digits = 4,
  file   = file.path(results_dir, "model_A37.doc")
)


###############################################################################
############################### Figure A1 #####################################
###############################################################################

#### Prepare colors
issues <- c(
  "Health","Labor","Education","Environment","Energy","Social welfare",
  "Macroeconomics","Law and crime","Domestic commerce","Defense","Foreign trade",
  "Civil rights","Agriculture","Immigration","Transport","Housing","Technology",
  "International affairs","Government operations","Public lands","EU","Culture"
)

colors <- c(
  rep("red", 6),
  rep("blue", 5),
  rep("grey", 11)
)

names(colors) <- issues

#### Plot
figure_a1 <- ggplot(data_figure_a1,
                    aes(x = CAP_issue, y = proportion, fill = CAP_issue)) +
  geom_bar(stat = "identity", position = "stack") +
  scale_fill_manual(values = colors) +
  facet_grid(. ~ party, scales = "free_x", space = "free_x") +
  theme_bw() +
  theme(axis.text.x = element_text(size = 18, angle = 90),
        axis.text.y = element_text(size = 18),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20),
        strip.text = element_text(size = 2),
        legend.position = "none") +
  labs(x = "\nPolicy issue", y = "Proportion\n")

#### Output
ggsave(
  filename = "figure_A1.pdf",
  plot = figure_a1,
  path = results_dir,
  width = 12,
  height = 8
)


###############################################################################
############################### Figure A2 #####################################
###############################################################################

#### Plot
figure_a2 <- ggplot(data_figure_a2,
                    aes(x = CAP_issue, y = proportion, fill = CAP_issue)) +
  geom_bar(stat = "identity", position = "stack") +
  scale_fill_manual(values = colors) +
  facet_grid(. ~ party, scales = "free_x", space = "free_x") +
  theme_bw() +
  theme(axis.text.x = element_text(size = 18, angle = 90),
        axis.text.y = element_text(size = 18),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20),
        strip.text = element_text(size = 20),
        legend.position = "none") +
  labs(x = "\nPolicy issue", y = "Share\n")

#### Output
ggsave(
  filename = "figure_A2.pdf",
  plot = figure_a2,
  path = results_dir,
  width = 12,
  height = 8
)


###############################################################################
############################### Figure A3 #####################################
###############################################################################

#### Plot
figure_a3 <- ggplot(data_figure_a3,
                    aes(x = engagement_count_on_initiation_day, y = n)) +
  geom_bar(stat = "identity") +
  geom_vline(xintercept = mean(data_figure_a3$mean_engagement),
             linetype = "dotted", color = "green", size = 2) +
  theme_minimal() +
  theme(axis.text.x = element_text(size = 20, angle = 90),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20),
        strip.text = element_text(size = 20),
        legend.position = "none") +
  labs(x = "\nEngagement ", y = "Frequency\n")

#### Output
ggsave(
  filename = "figure_A3.pdf",
  plot = figure_a3,
  path = results_dir,
  width = 10,
  height = 8
)


###############################################################################
############################### Figure A4 #####################################
###############################################################################

#### Plot
figure_a4 <- ggplot(data_figure_a4,
                    aes(x = news_articles_day_after, y = n)) +
  geom_bar(stat = "identity") +
  geom_vline(xintercept = mean(data_figure_a4$mean_news_articles),
             linetype = "dotted", color = "green", size = 2) +
  theme_minimal() +
  theme(axis.text.x = element_text(size = 20, angle = 90),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20),
        strip.text = element_text(size = 20),
        legend.position = "none") +
  labs(x = "\nNews media agenda", y = "Frequency\n")

#### Output
ggsave(
  filename = "figure_A4.pdf",
  plot = figure_a4,
  path = results_dir,
  width = 10,
  height = 8
)


###############################################################################
############################### Figure A5 #####################################
###############################################################################

#### Plot
figure_a5 <- data_figure_a5 %>%
  ggplot(aes(x = period, y = irf_elevation_reaction,
             ymin = lower_elevation_reaction, ymax = upper_elevation_reaction)) +
  geom_hline(yintercept = 0, color = "red") +
  geom_ribbon(fill = "grey", alpha = 0.2) +
  geom_line() +
  theme_light() +
  ylab("Number of tweets by competing party MPs\n") +
  xlab("\nNumber of minutes since the one-tweet increase in elevation") +
  theme_minimal() +
  theme(axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 22),
        axis.title.y = element_text(size = 22),
        strip.text = element_text(size = 22)) +
  scale_x_continuous(breaks = data_figure_a5$period)

#### Output
ggsave(
  filename = "figure_A5.pdf",
  plot = figure_a5,
  path = results_dir,
  width = 12,
  height = 8
)


###############################################################################
############################### Figure A6 #####################################
###############################################################################

#### Plot
figure_a6 <- data_figure_a6 %>%
  ggplot(aes(x = period, y = irf_elevation_reaction,
             ymin = lower_elevation_reaction, ymax = upper_elevation_reaction)) +
  geom_hline(yintercept = 0, color = "red") +
  geom_ribbon(fill = "grey", alpha = 0.2) +
  geom_line() +
  theme_light() +
  ylab("Cumulated number of tweets by competing party MPs\n") +
  xlab("\nNumber of minutes since the one-tweet increase in elevation") +
  theme_minimal() +
  theme(axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 18),
        axis.title.y = element_text(size = 18),
        strip.text = element_text(size = 22)) +
  scale_x_continuous(breaks = data_figure_a6$period)

#### Output
ggsave(
  filename = "figure_A6.pdf",
  plot = figure_a6,
  path = results_dir,
  width = 12,
  height = 8
)


###############################################################################
############################### Figure A7 #####################################
###############################################################################

#### Plot
figure_a7 <- data_figure_a7 %>%
  ggplot(aes(x = period, y = irf_elevation_reaction,
             ymin = lower_elevation_reaction, ymax = upper_elevation_reaction)) +
  geom_hline(yintercept = 0, color = "red") +
  geom_ribbon(fill = "grey", alpha = 0.2) +
  geom_line() +
  theme_light() +
  ylab("Number of tweets by competing party MPs\n") +
  xlab("\nNumber of minutes since the one-tweet increase in elevation") +
  theme_minimal() +
  theme(axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 22),
        axis.title.y = element_text(size = 22),
        strip.text = element_text(size = 22)) +
  scale_x_continuous(breaks = data_figure_a7$period)

#### Output
ggsave(
  filename = "figure_A7.pdf",
  plot = figure_a7,
  path = results_dir,
  width = 12,
  height = 8
)


###############################################################################
############################### Figure A8 #####################################
###############################################################################

#### Plot
figure_a8 <- data_figure_a8 %>%
  ggplot(aes(x = period, y = irf_elevation_reaction,
             ymin = lower_elevation_reaction, ymax = upper_elevation_reaction)) +
  geom_hline(yintercept = 0, color = "red") +
  geom_ribbon(fill = "grey", alpha = 0.2) +
  geom_line() +
  theme_light() +
  ylab("Cumulated number of tweets by competing party MPs\n") +
  xlab("\nNumber of minutes since the one-tweet increase in elevation") +
  theme_minimal() +
  theme(axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 18),
        axis.title.y = element_text(size = 18),
        strip.text = element_text(size = 22)) +
  scale_x_continuous(breaks = data_figure_a8$period)

#### Output
ggsave(
  filename = "figure_A8.pdf",
  plot = figure_a8,
  path = results_dir,
  width = 12,
  height = 8
)


###############################################################################
############################### Figure A9 #####################################
###############################################################################

#### Plot
figure_a9 <- ggplot(data_figure_a9,
                    aes(x = who, y = count, fill = who)) +
  geom_bar(stat = "identity") +
  labs(x = "\nWho created most tweets",
       y = "Number of initiation days\n") +
  theme(axis.text.x = element_text(size = 20, angle = 25, hjust = 1),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 22),
        axis.title.y = element_text(size = 22),
        strip.text   = element_text(size = 22),
        plot.title   = element_text(size = 22),
        legend.position = "none",
        legend.text  = element_text(size = 20)) +
  geom_text(aes(label = count), vjust = 0.4, size = 5.5) +
  facet_wrap(~ hour, ncol = 2)

#### Output
ggsave(
  filename = "figure_A9.pdf",
  plot = figure_a9,
  path = results_dir,
  width = 14,
  height = 8
)


###############################################################################
############################### Figure A10 ####################################
###############################################################################

#### Plot
figure_a10 <- ggplot(data_figure_a10,
                     aes(x = who, y = count, fill = who)) +
  geom_bar(stat = "identity") +
  labs(x = "\nWho created most tweets",
       y = "Number of initiation days\n") +
  theme(axis.text.x = element_text(size = 16, angle = 25, hjust = 1),
        axis.text.y = element_text(size = 16),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20),
        strip.text   = element_text(size = 20),
        plot.title   = element_text(size = 20),
        legend.position = "none",
        legend.text  = element_text(size = 16)) +
  geom_text(aes(label = count), vjust = 0.4, size = 5.5) +
  facet_wrap(~ hour, ncol = 2)

#### Output
ggsave(
  filename = "figure_A10.pdf",
  plot = figure_a10,
  path = results_dir,
  width = 12,
  height = 8
)


###############################################################################
############################### Figure A11 ####################################
###############################################################################

#### Plot
figure_a11 <- ggplot(data_figure_a11,
                     aes(x = who, y = average, fill = who)) +
  geom_bar(stat = "identity") +
  labs(x = "\nActors",
       y = "Number of tweets at the specific time\n") +
  theme(axis.text.x = element_text(size = 16, angle = 25, hjust = 1),
        axis.text.y = element_text(size = 16),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20),
        strip.text   = element_text(size = 20),
        plot.title   = element_text(size = 20),
        legend.position = "none",
        legend.text  = element_text(size = 20)) +
  geom_text(aes(label = scales::number(average, accuracy = 0.01)),
            vjust = 0.4, size = 5.5) +
  facet_wrap(~ hour, ncol = 2)

#### Output
ggsave(
  filename = "figure_A11.pdf",
  plot = figure_a11,
  path = results_dir,
  width = 12,
  height = 8
)


###############################################################################
############################### Figure A12 ####################################
###############################################################################

#### Plot
figure_a12 <- ggplot(data_figure_a12,
                     aes(x = who, y = average, fill = who)) +
  geom_bar(stat = "identity") +
  labs(x = "\nActors",
       y = "Number of posts at the specific time\n") +
  theme(axis.text.x = element_text(size = 16, angle = 25, hjust = 1),
        axis.text.y = element_text(size = 16),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20),
        strip.text   = element_text(size = 20),
        plot.title   = element_text(size = 20),
        legend.position = "none",
        legend.text  = element_text(size = 16)) +
  geom_text(aes(label = scales::number(average, accuracy = 0.01)),
            vjust = 0.4, size = 5.5) +
  facet_wrap(~ hour, ncol = 2)

#### Output
ggsave(
  filename = "figure_A12.pdf",
  plot = figure_a12,
  path = results_dir,
  width = 12,
  height = 8
)


###############################################################################
############################### Figure A13 ####################################
###############################################################################

#### Plot
figure_a13 <- data_figure_a13 %>%
  ggplot(aes(x = period, y = irf_elevation_reaction,
             ymin = lower_elevation_reaction, ymax = upper_elevation_reaction)) +
  geom_hline(yintercept = 0, color = "red") +
  geom_ribbon(fill = "grey", alpha = 0.2) +
  geom_line() +
  theme_light() +
  ylab("Number of tweets by competing party MPs\n") +
  xlab("\nNumber of minutes since the one-tweet increase in elevation") +
  theme_minimal() +
  theme(axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 22),
        axis.title.y = element_text(size = 22),
        strip.text   = element_text(size = 22)) +
  scale_x_continuous(breaks = data_figure_a13$period)

#### Output
ggsave(
  filename = "figure_A13.pdf",
  plot = figure_a13,
  path = results_dir,
  width = 12,
  height = 8
)


###############################################################################
############################### Figure A14 ####################################
###############################################################################

#### Plot
figure_a14 <- data_figure_a14 %>%
  ggplot(aes(x = period, y = irf_elevation_reaction,
             ymin = lower_elevation_reaction, ymax = upper_elevation_reaction)) +
  geom_hline(yintercept = 0, color = "red") +
  geom_ribbon(fill = "grey", alpha = 0.2) +
  geom_line() +
  theme_light() +
  ylab("Number of tweets by competing party MPs\n") +
  xlab("\nNumber of minutes since the one-tweet increase in elevation") +
  theme_minimal() +
  theme(axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 22),
        axis.title.y = element_text(size = 22),
        strip.text   = element_text(size = 22)) +
  scale_x_continuous(breaks = data_figure_a14$period)

#### Output
ggsave(
  filename = "figure_A14.pdf",
  plot = figure_a14,
  path = results_dir,
  width = 12,
  height = 8
)


###############################################################################
############################### Figure A15 ####################################
###############################################################################

#### Plot
figure_a15 <- ggplot(data_figure_a15,
                     aes(x = mp_gender, y = mean_elevation_rate)) +
  geom_col(fill = "black", alpha = 0.7) +
  facet_wrap(~ country, scales = "free_y") +
  labs(x = "Gender",
       y = "Mean Elevation Rate (%)") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 12),
        axis.text.y = element_text(size = 12),
        axis.title  = element_text(size = 14),
        strip.text  = element_text(size = 16),
        plot.title  = element_text(size = 16))

#### Output
ggsave(
  filename = "figure_A15.pdf",
  plot = figure_a15,
  path = results_dir,
  width = 12,
  height = 8
)


###############################################################################
############################### Figure A16 ####################################
###############################################################################

#### Plot
figure_a16 <- ggplot(data_figure_a16,
                     aes(x = seniority_group, y = mean_elevation_rate)) +
  geom_col(fill = "black", alpha = 0.7) +
  facet_wrap(~ country, scales = "free_y") +
  labs(x = "Parliamentary Seniority",
       y = "Mean Elevation Rate (%)") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 12),
        axis.text.y = element_text(size = 12),
        axis.title  = element_text(size = 14),
        strip.text  = element_text(size = 16),
        plot.title  = element_text(size = 16))

#### Output
ggsave(
  filename = "figure_A16.pdf",
  plot = figure_a16,
  path = results_dir,
  width = 12,
  height = 8
)


###############################################################################
############################### Figure A17 ####################################
###############################################################################

#### Plot
figure_a17 <- ggplot(data_figure_a17,
                     aes(x = minister_status, y = mean_elevation_rate)) +
  geom_col(fill = "black", alpha = 0.7) +
  facet_wrap(~ country, scales = "free_y") +
  labs(x = "Ministerial Status",
       y = "Mean Elevation Rate (%)") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 12),
        axis.text.y = element_text(size = 12),
        axis.title  = element_text(size = 14),
        strip.text  = element_text(size = 16),
        plot.title  = element_text(size = 16))

#### Output
ggsave(
  filename = "figure_A17.pdf",
  plot = figure_a17,
  path = results_dir,
  width = 12,
  height = 8
)


###############################################################################
############################### Figure A18 ####################################
###############################################################################

#### Plot
figure_a18 <- data_figure_a18 %>%
  dplyr::mutate(country = dplyr::case_when(
    country == "dk" ~ "Denmark",
    country == "uk" ~ "United Kingdom",
    TRUE ~ country
  )) %>%
  dplyr::group_by(country, mp_gender) %>%
  dplyr::summarise(mean_tweets = mean(tweet_count), .groups = "drop") %>%
  ggplot(aes(x = mp_gender, y = mean_tweets)) +
  geom_col(fill = "grey30") +
  facet_wrap(~ country, scales = "free") +
  labs(x = "Gender",
       y = "Mean Tweet Count") +
  theme_minimal() +
  theme(axis.text.x = element_text(size = 14),
        axis.text.y = element_text(size = 12),
        axis.title  = element_text(size = 14),
        strip.text  = element_text(size = 16),
        plot.title  = element_text(size = 16)) +
  theme(legend.position = "none")

#### Output
ggsave(
  filename = "figure_A18.pdf",
  plot = figure_a18,
  path = results_dir,
  width = 12,
  height = 8
)


###############################################################################
############################### Figure A19 ####################################
###############################################################################

#### Plot
figure_a19 <- data_figure_a19 %>%
  dplyr::mutate(country = dplyr::case_when(
    country == "dk" ~ "Denmark",
    country == "uk" ~ "United Kingdom",
    TRUE ~ country
  )) %>%
  dplyr::group_by(country, seniority_group) %>%
  dplyr::summarise(
    mean_tweets = mean(tweet_count, na.rm = TRUE),
    n = dplyr::n(),
    .groups = "drop"
  ) %>%
  ggplot(aes(x = seniority_group, y = mean_tweets)) +
  geom_col(fill = "black", alpha = 0.7) +
  facet_wrap(~ country, scales = "free_y") +
  labs(x = "Parliamentary Seniority",
       y = "Mean Elevation Count") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 12),
        axis.text.y = element_text(size = 12),
        axis.title  = element_text(size = 14),
        strip.text  = element_text(size = 16),
        plot.title  = element_text(size = 16))

#### Output
ggsave(
  filename = "figure_A19.pdf",
  plot = figure_a19,
  path = results_dir,
  width = 12,
  height = 8
)


###############################################################################
############################### Figure A20 ####################################
###############################################################################

#### Plot
figure_a20 <- data_figure_a20 %>%
  dplyr::group_by(country, minister_status) %>%
  dplyr::summarise(
    mean_tweets = mean(tweet_count, na.rm = TRUE),
    n = dplyr::n(),
    .groups = "drop"
  ) %>%
  ggplot(aes(x = minister_status, y = mean_tweets)) +
  geom_col(fill = "black", alpha = 0.7) +
  facet_wrap(~ country, scales = "free_y") +
  labs(x = "Ministerial Status",
       y = "Mean Elevation Count") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 12),
        axis.text.y = element_text(size = 12),
        axis.title  = element_text(size = 14),
        strip.text  = element_text(size = 16),
        plot.title  = element_text(size = 16))

#### Output
ggsave(
  filename = "figure_A20.pdf",
  plot = figure_a20,
  path = results_dir,
  width = 12,
  height = 8
)


###############################################################################
############################### Figure A21 ####################################
###############################################################################

#### Plot
figure_a21 <- ggplot(data_figure_a21,
                     aes(x = date, y = count, fill = party)) +
  geom_bar(stat = "identity") +
  scale_fill_manual(values = c("Enhedslisten" = "red",
                               "Dansk Folkeparti" = "blue"),
                    na.value = "gray") +
  theme_minimal() +
  theme(axis.text.x = element_text(size = 26, angle = 90),
        axis.text.y = element_text(size = 28),
        axis.title.x = element_text(size = 28),
        axis.title.y = element_text(size = 28),
        strip.text   = element_text(size = 28),
        legend.position = "none") +
  labs(x = "\nDate", y = "Frequency\n") +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 5))

ggsave(
  filename = "figure_A21.pdf",
  plot = figure_a21,
  path = results_dir,
  width = 30,
  height = 12,
  units = "in",
  dpi = 300
)


###############################################################################
############################### Figure A22 ####################################
###############################################################################

#### Plot
figure_a22 <- ggplot(
  data_figure_a22,
  aes(x = days_before_tweet,
      y = average_number_tweets_other_parties,
      group = tweet_type,
      linetype = tweet_type)
) +
  geom_line(size = 1.5) +
  geom_point(size = 3) +
  labs(
    x = "\nDays before the post",
    y = "Average number of tweets by MPs from all other parties\n"
  ) +
  theme_minimal() +
  ylim(0, 20) +
  scale_color_manual(values = "black") +
  scale_linetype_manual(values = c("Initiations" = "solid",
                                   "Non-Initiations" = "dashed")) +
  theme(
    plot.title = element_text(size = 18, hjust = 0.5),
    text       = element_text(size = 18),
    legend.position = "none",
    axis.text.x = element_text(size = 15, margin = margin(t = 20), angle = 25),
    axis.text.y = element_text(size = 15, margin = margin(r = 20)),
    axis.title  = element_text(size = 15, margin = margin(b = 20))
  )

#### Output
ggsave(
  filename = "figure_A22.pdf",
  plot = figure_a22,
  path = results_dir,
  width = 12,
  height = 8
)


###############################################################################
############################### Figure A23 ####################################
###############################################################################

#### Plot
figure_a23 <- ggplot(
  data_figure_a23,
  aes(x = days_before_tweet,
      y = average_number_tweets_other_parties,
      group = tweet_type,
      linetype = tweet_type)
) +
  geom_line(size = 1.5) +
  geom_point(size = 3) +
  labs(
    x = "\nDays before the post",
    y = "Average number of tweets by MPs from all other parties\n"
  ) +
  theme_minimal() +
  ylim(0, 20) +
  scale_color_manual(values = "black") +
  scale_linetype_manual(values = c("Initiations" = "solid",
                                   "Non-Initiations" = "dashed")) +
  theme(
    plot.title = element_text(size = 18, hjust = 0.5),
    text       = element_text(size = 18),
    legend.position = "none",
    axis.text.x = element_text(size = 15, margin = margin(t = 20), angle = 25),
    axis.text.y = element_text(size = 15, margin = margin(r = 20)),
    axis.title  = element_text(size = 15, margin = margin(b = 20))
  )

#### Output
ggsave(
  filename = "figure_A23.pdf",
  plot = figure_a23,
  path = results_dir,
  width = 12,
  height = 8
)


###############################################################################
############################### Figure A24 ####################################
###############################################################################

#### Plot
figure_a24 <- ggplot(
  data_figure_a24,
  aes(x = days_before_tweet,
      y = average_news_articles,
      group = tweet_type,
      linetype = tweet_type)
) +
  geom_line(size = 1.5) +
  geom_point(size = 3) +
  labs(
    x = "\nDays before the post",
    y = "Average number of news media articles\n"
  ) +
  theme_minimal() +
  ylim(0, 10) +
  scale_color_manual(values = "black") +
  scale_linetype_manual(values = c("Initiations" = "solid",
                                   "Non-Initiations" = "dashed")) +
  theme(
    plot.title = element_text(size = 18, hjust = 0.5),
    text       = element_text(size = 18),
    legend.position = "none",
    axis.text.x = element_text(size = 18, margin = margin(t = 20), angle = 25),
    axis.text.y = element_text(size = 18, margin = margin(r = 20)),
    axis.title  = element_text(size = 18, margin = margin(b = 20))
  )

#### Output
ggsave(
  filename = "figure_A24.pdf",
  plot = figure_a24,
  path = results_dir,
  width = 12,
  height = 8
)


###############################################################################
############################### Figure A25 ####################################
###############################################################################

#### Plot
figure_a25 <- ggplot(
  data_figure_a25,
  aes(x = day, y = average_elevation_tweets, group = tweet_type)
) +
  geom_line(aes(linetype = tweet_type), size = 1.0) +
  geom_point(aes(color = tweet_type), size = 3) +
  labs(
    x = "Day",
    y = "Average number of tweets about the issue of the post",
    title = ""
  ) +
  theme_minimal() +
  scale_color_manual(values = c("Non-initiation" = "black",
                                "Initiation" = "black")) +
  scale_linetype_manual(values = c("Non-initiation" = "dashed",
                                   "Initiation" = "solid")) +
  scale_x_discrete(limits = rev(levels(data_figure_a25$day))) +
  theme(
    plot.title = element_text(size = 20, hjust = 0.5),
    text       = element_text(size = 18),
    axis.text.x = element_text(size = 16, margin = margin(t = 20)),
    axis.text.y = element_text(size = 16, margin = margin(r = 20)),
    axis.title  = element_text(size = 18, margin = margin(b = 20)),
    legend.position = "none"
  ) +
  geom_text(
    aes(x = day, label = round(average_elevation_tweets, 1)),
    vjust = -0.5,
    size = 6
  )

#### Output
ggsave(
  filename = "figure_A25.pdf",
  plot = figure_a25,
  path = results_dir,
  width = 12,
  height = 8
)


###############################################################################
############################### Figure A26 ####################################
###############################################################################

#### Plot
figure_a26 <- data_figure_a26 %>%
  ggplot(aes(x = period, y = irf_elevation_reaction,
             ymin = lower_elevation_reaction, ymax = upper_elevation_reaction)) +
  geom_hline(yintercept = 0, color = "red") +
  geom_ribbon(fill = "grey", alpha = 0.2) +
  geom_line() +
  theme_light() +
  ylab("Number of tweets by competing party MPs\n") +
  xlab("\nNumber of minutes since the one-tweet increase in elevation") +
  theme_minimal() +
  theme(axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20),
        strip.text   = element_text(size = 20)) +
  scale_x_continuous(breaks = data_figure_a26$period)

#### Output
ggsave(
  filename = "figure_A26.pdf",
  plot = figure_a26,
  path = results_dir,
  width = 12,
  height = 8
)


###############################################################################
############################### Table A2 ######################################
###############################################################################

#### Output
write.csv(
  data_table_a2,
  file = file.path(results_dir, "Table_A2.csv"),
  row.names = FALSE
)
